If there are errors take place, you can run install.packages({missing package name}) to install packages.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ ggplot2 3.2.1 ✔ purrr 0.3.3
## ✔ tibble 3.0.3 ✔ dplyr 1.0.2
## ✔ tidyr 1.1.2 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:purrr':
##
## transpose
library(sp)
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(table1)
##
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
##
## units, units<-
library(tm)
## Loading required package: NLP
##
## Attaching package: 'NLP'
## The following object is masked from 'package:ggplot2':
##
## annotate
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## Registered S3 method overwritten by 'spatstat':
## method from
## print.boxx cli
##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.64-1 is out of date by more than 8 months; we recommend upgrading to the latest version.
##
## Attaching package: 'spatstat'
## The following object is masked from 'package:data.table':
##
## shift
library(here)
## here() starts at /Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment
library(sp)
library(rgeos)
## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.7.2-CAPI-1.11.2
## Linking to sp version: 1.4-2
## Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
##
## Attaching package: 'maptools'
## The following object is masked from 'package:table1':
##
## label
library(tmap)
library(sf)
library(geojson)
##
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
##
## polygon
library(geojsonio)
##
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
##
## pretty
library(tmaptools)
library(RColorBrewer)
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
The path below is mine, you should set your own work path
setwd("/Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment/")
What you should keep in mind is that this shape file should be run in the complete ESRI dir
# London_Borough=st_read("F:/spat/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")
London_Borough = st_read("dataset/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")
## Reading layer `London_Borough_Excluding_MHW' from data source `/Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment/Dataset/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## CRS: 27700
# plot the map
plot(st_geometry(London_Borough))
df = fread("Dataset/London_NCR_GSS_Added.csv")
#3.1 select the years
df$year = year(df$dateCreated)
table(df$year)
##
## 2017 2018 2019 2020 2021
## 1 4 696 362 3
#select 2019 & 2020
df = df[year==2019|year==2020,]
table(df$year)
##
## 2019 2020
## 696 362
#3.2 show the table of the distribution of two years samples
table1(~county|factor(year),data=df)
| 2019 (N=696) |
2020 (N=362) |
Overall (N=1058) |
|
|---|---|---|---|
| county | |||
| London | 150 (21.6%) | 150 (41.4%) | 300 (28.4%) |
| London Borough of Camden | 87 (12.5%) | 179 (49.4%) | 266 (25.1%) |
| London Borough of Ealing | 65 (9.3%) | 1 (0.3%) | 66 (6.2%) |
| London Borough of Greenwich | 2 (0.3%) | 0 (0%) | 2 (0.2%) |
| London Borough of Hackney | 62 (8.9%) | 0 (0%) | 62 (5.9%) |
| London Borough of Hammersmith and Fulham | 2 (0.3%) | 3 (0.8%) | 5 (0.5%) |
| London Borough of Hounslow | 78 (11.2%) | 0 (0%) | 78 (7.4%) |
| London Borough of Islington | 3 (0.4%) | 15 (4.1%) | 18 (1.7%) |
| London Borough of Lambeth | 118 (17.0%) | 2 (0.6%) | 120 (11.3%) |
| London Borough of Richmond upon Thames | 8 (1.1%) | 9 (2.5%) | 17 (1.6%) |
| London Borough of Southwark | 1 (0.1%) | 0 (0%) | 1 (0.1%) |
| London Borough Of Southwark | 57 (8.2%) | 0 (0%) | 57 (5.4%) |
| London Borough of Waltham Forest | 60 (8.6%) | 2 (0.6%) | 62 (5.9%) |
| London Borough of Wandsworth | 3 (0.4%) | 1 (0.3%) | 4 (0.4%) |
head(df)
## chargeDeviceID latitude longitude town county
## 1: b353888fbc3ee702b41d30669a23e12d 51.54491 -0.00876 Stratford London
## 2: f53d6d41d4ea910e11d4ea914d58b803 51.54491 -0.00876 Stratford London
## 3: 61ef7c5a4cdc44aa15d46293f1f185b8 51.54491 -0.00876 Stratford London
## 4: 3ae50acb1eb21a62b4f1d6620f155c10 51.54491 -0.00876 Stratford London
## 5: 3c348a8f6ea52056921b931828359346 51.54491 -0.00876 Stratford London
## 6: 6e785fdfb123ff8b192f26cf595f74b1 51.54491 -0.00876 Stratford London
## postcode chargeDeviceStatus dateCreated dateUpdated
## 1: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 2: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 3: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 4: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 5: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 6: E20 1YY In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## lastUpdated locationType GSS_CODE year
## 1: Public car park E09000025 2020
## 2: Public car park E09000025 2020
## 3: Public car park E09000025 2020
## 4: Public car park E09000025 2020
## 5: Public car park E09000025 2020
## 6: Public car park E09000025 2020
#5.1 process the data to divide them by years(2019 & 2020)
df1<-subset(df,year==2019) #2019 year
df2<-subset(df,year==2020) #2020 year
# EV charge points created in 2019
# sdf1<-merge(London_Borough,df1,by="GSS_CODE")
sdf1<-merge(London_Borough,df1,by="GSS_CODE",all = TRUE)
sdf1<-sdf1[,c("GSS_CODE","geometry","longitude","latitude")]
# EV charge points created in 2020
# sdf2<-merge(London_Borough,df2,by="GSS_CODE")
sdf2<-merge(London_Borough,df2,by="GSS_CODE",all = TRUE)
sdf2<-sdf2[,c("GSS_CODE","geometry","longitude","latitude")]
### 4.2 The density of data in 2019
# Data preparation
nsdf1 = sdf1%>%
add_count(GSS_CODE)%>%
mutate(area=st_area(.))%>%
# Use dplyr::mutate to calculate the density of the charge point for each borough
mutate(density=n*1000*1000/area)
# because the st_area default unit is square metre
# select the following variables---"density","GSS_CODE","n"(the count of GSS_CODE)
nsdf1 = dplyr::select(nsdf1,density,GSS_CODE, n)
nsdf1 = nsdf1%>%
group_by(GSS_CODE)%>%
summarise(density =first(density),GSS_CODE=first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)
tmap_mode("plot")
## tmap mode set to plotting
# plot the figure2: The distribution of the density of the London charge points in 2019
tm_shape( nsdf1) +
tm_compass( north = 0,
type = "4star",
text.size = 0.8,
size = 2.5,
show.labels = 1,
cardinal.directions = c("N", "E", "S", "W"),
lwd = 1,
position = c("left","top"),
bg.color = NA,
bg.alpha = NA,
just = NA,
fontsize = 1.5) +
tm_polygons("density",
style="jenks",
palette="RdPu",
midpoint=NA,
popup.vars=c("GSS_CODE", "density"),
title="Density per square kilometre (2019)"
)
## Warning: The argument fontsize of tm_compass is deprecated. It has been renamed
## to text.size
Figure the London boroughs map with their names
### 4.3 The density of data in 2020
nsdf2 = sdf2%>%
add_count(GSS_CODE)%>%
mutate(area=st_area(.))%>%
mutate(units::set_units(area,km^2))%>%
#then density of the points per ward
mutate(density=n*1000*1000/area)
nsdf2 = dplyr::select(nsdf2,density,GSS_CODE, n)
nsdf2 = nsdf2%>%
group_by(GSS_CODE)%>%
summarise(density =first(density),GSS_CODE=first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)
tmap_mode("plot")
## tmap mode set to plotting
tm_shape( nsdf2) +
tm_compass( north = 0,
type = "4star",
text.size = 0.8,
size = 2.5,
show.labels = 1,
cardinal.directions = c("N", "E", "S", "W"),
lwd = 1,
position = c("left","top"),
bg.color = NA,
bg.alpha = NA,
just = NA,
fontsize = 1.5) +
tm_polygons("density",
style="jenks",
palette="PuOr",
midpoint=NA,
popup.vars=c("GSS_CODE", "density"),
title="Density per square kilometre (2020)")
## Warning: The argument fontsize of tm_compass is deprecated. It has been renamed
## to text.size
Since the sample in 2019 & 2020 is too small to run a good result, I analysis the autocorrelation based on the samples which contain these two yearsl
###5.1 generate the data for analysis
sdf = merge(London_Borough,df,by="GSS_CODE",all = TRUE)
sdf = sdf[,c("GSS_CODE","geometry","longitude","latitude")]
nsdf = sdf%>%
add_count(GSS_CODE)%>%
mutate(area=st_area(.))%>%
#then density of the points per ward
mutate(density=n*1000*1000/area)
nsdf = dplyr::select(nsdf,density,GSS_CODE, n)
nsdf = nsdf%>%
group_by(GSS_CODE)%>%
summarise(density = first(density), GSS_CODE = first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)
###5.2 First plot the centroids of all boroughs in London
coordsW = nsdf%>%
st_centroid()%>%
st_geometry()
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
plot(coordsW,axes=TRUE)
LWard_nb <- nsdf%>%poly2nb(.,queen=T)
plot the neighbours list we creat
plot(nsdf$geometry)
plot(LWard_nb, st_geometry(coordsW), col="blue", add = T)
# plot(LWard_nb, st_geometry(coordsW), col="red")
# plot(nsdf$geometry, add=T)
# plot(LWard_nb, st_geometry(coordsW), col="red")
# jpeg(file="saving_plot1.jpeg")
# dev.off()
### 5.3 create a spatial weights object from these weights
Lward.lw = nb2listw(LWard_nb, style="C")
head(Lward.lw$neighbours)
## [[1]]
## [1] 7 12 19 30 33
##
## [[2]]
## [1] 16 25 26
##
## [[3]]
## [1] 5 7 10 14 15
##
## [[4]]
## [1] 6 11
##
## [[5]]
## [1] 3 7 9 13 15 20 33
##
## [[6]]
## [1] 4 8 11 22 23 28
### 5.4 Calculate the Global Moran’I Index
I_LWard_Global_Density = nsdf %>%
pull(density) %>%
as.vector()%>%
moran.test(.,Lward.lw)
names(I_LWard_Global_Density)
## [1] "statistic" "p.value" "estimate" "alternative" "method"
## [6] "data.name"
head(I_LWard_Global_Density)
## $statistic
## Moran I statistic standard deviate
## 0.03513093
##
## $p.value
## [1] 0.4859877
##
## $estimate
## Moran I statistic Expectation Variance
## -0.028283348 -0.031250000 0.007131055
##
## $alternative
## [1] "greater"
##
## $method
## [1] "Moran I test under randomisation"
##
## $data.name
## [1] ". \nweights: Lward.lw \n"
I_LWard_Local_Density = nsdf %>%
pull(density) %>%
as.vector()%>%
localmoran(., Lward.lw)%>%
as_tibble()
nsdf<-nsdf%>%
mutate(density_I = as.numeric(I_LWard_Local_Density$Ii))%>%
mutate(density_Iz =as.numeric(I_LWard_Local_Density$Z.Ii))
tmap_mode("view")
## tmap mode set to interactive viewing
#set the group and colour
summary(nsdf$density_Iz)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.813183 -0.007425 0.133178 0.046786 0.449411 0.975067
breaks1 = seq(-3,1,0.5)
# breaks2 = c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
# Depends on the max and min value in Moran's I
MoranColours = rev(brewer.pal(8, "RdGy"))
# Plot the map
tm_shape(nsdf) +
tm_polygons("density_Iz",
style="fixed",
breaks=breaks1,
palette=MoranColours,
midpoint=NA,
title="Local Moran's I,EV charge points in London")
### 5.6 GI score
Gi_LWard_Local_Density = nsdf %>%
pull(density) %>%
as.vector()%>%
localG(., Lward.lw)
head(Gi_LWard_Local_Density)
## [1] 2.3573180 -0.7667792 1.8494548 -0.6694542 1.1946019 0.6283400
Gi_nsdf = nsdf %>%
mutate(density_G = as.numeric(Gi_LWard_Local_Density))
GIColours = rev(brewer.pal(8, "RdBu"))
#now plot on an interactive map
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(Gi_nsdf) +
tm_polygons("density_G",
style="fixed",
breaks=breaks1,
palette=GIColours,
midpoint=NA,
title="Gi*, EV charge points in London")
## Warning: Values have found that are higher than the highest break
Thanks for watching!